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Abstract. In this paper we investigate how the inclusion of scattering of the stellar radiation into a passive flaring 
disk model affects its structure and spectral energy distribution, and whether neglecting it could significantly 
decrease the model reliability. In order to address these questions we construct a detailed 1+1D vertical structure 
model in which the scattering properties of the dust can be varied. Models are presented with and without dust 
scattering, and for different albedos and phase functions. It is found that scattering has the effect of reducing 
the disk temperature at all heights, so that the disk "shrinks", i.e., the the density at all intermediate heights 
decreases. However, this effect in most cases is more than compensated by the increase of the total extinction 
(absorption + scattering) cross section, so that the surface scale height increases, and images in scattered light will 
see a slightly thicker disk. The integrated infrared emission decreases as the albedo increases, because an increasing 
part of the flux captured by the disk is reflected away instead of absorbed and reprocessed. The reduction of the 
infrared thermal emission of the disk is stronger at short wavelengths (near infrared) and practically negligible at 
millimeter wavelengths. For relatively low albedo (u < 0.5), or for strongly forward-peaked scattering (g roughly 
> 0.8), the infrared flux reduction is relatively small. 
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1. Introduction 

In recent years, models of irradiated passive protoplane- 
tary disks have proven to be quite successful in explaining 
the spectral energy distributions (SEDs) of T Tauri (TTS) 
and H erbig Ae/Be (HAeBe) st ars (Kenyon & Hartmann 
119871 Calvet et al. Jl99lll992t Malbet & B ertout 
Chiang et al. Il997l200lt D'Alessi o et al. Il99. " 
Malbet et 



1991 



1999: 



al. 1200 It Natta et al. l2001t Dullemond et 
al. l200lt Dominik et al. l2003j) . Though the exact geometry 
of such disks is still an issue of debate, the basic idea of a 
passively reprocessing dusty disk, possibly with a remnant 
tenuous envelope surrounding it, thusfar seems to with- 
stand the test of time. However, it is clear that analyzing 
the SEDs alone is in most cases not sufficient to constrain 
th e disk 's parameters or to test theories (e.g. Chiang et 
al. l200llh and that spatially resolved images at various 
wavelengths are required in addition. Spatially-resolved 
images at millimeter wavelengths and optical and near- 
infrared images in scattered light provide valuable con- 
straints to the disk geometry and to the properties of the 
dust in thes e disks (see, for example, D'Alessio et al. l200lt 
Testi et al. l200ll 120031 for the use of millimeter maps to 



constrain outer disk properties and Cotera et al. l200ll and 
references therein for the interpretation of disk images 
in scattered light). Models that self-consistently and si- 
multaneously compute the disk structure, the SED and 
scattered-light images are required in order to fully ex- 
ploit all the observed data. Some self-c onsist ent models 
of this k ind exist (e.g. D'Alessio et al. l200lt Malbet et 
al. l200l[h but many other models use the assumption of 
zero albedo while computing the structure and SED of 
the disk (e.g. the models of Chiang et al. and those of 
Dullemond et al.). Most of the calculations of scattered 
light disk images use parametric disk models to de scribe 
the d isk structu re (se e, for example, Wood et al. [1.998, 
l2002t Wolf et al. l2003t) . Although very successful in deriv- 
ing a first characterization of the disk surface properties, 
this approach does not take into account the effects of 
scattering on the structure and thermal emission of the 
disk. 

In this paper, we explore systematically the effect on 
the disk structure and emission of including scattering 
of the stellar radiation by the grains in the disk. We do 
this by varying over a large range of values the albedo 
and phase function of the grains, and we consider disks 
heated by a typical Herbig Ae stars and by a typical TTS. 
Ou r stud y is somewhat complementary to D'Alessio et 
al. l|200ljh who present disk models where the grain size 
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distribution is extended to very large grains and the cross 
sections for absorption and scattering are changed accord- 
ingly. Their approach is more physical, since it relates the 
adopted grain model to the specific physical process of 
grain growth. However, a parametric study like ours has 
the advantage of greater freedom in choosing the dust 
properties, and this is particularly important when one 
is interested in understanding the relevance of scattering 
in disk models in general. 

In §2 we present a self-consistent 1+1D vertical struc- 
ture model for an irradiated passive flaring disk around a 
TTS or a HAeBe star in which the scattering properties 
of the dust grains can be varied freely. In the case of zero 
albedo this model reduces to the earlier publi shed model 
by Dullemond, van Zadelhoff & Natta (|2002t henceforth 
DZN02). For non-zero albedo the differences to the latter 
model can be investigated in detail, and an assessment can 
be made of the role that scattering of the stellar radiation 
plays in the structure and SED of these disks. We illus- 
trate the results for a disk annulus in §3, and for full disk 
models in §4. A summary and conclusions follow in §5. 



2. The vertical structure model with scattering 

The disk vertical structure model that we present in this 
section is based for most part on the 1+1D model of 
DZN02. The continuum radiative transfer, which is re- 
quired in order to compute the dust temperature, is done 
in two stages. In stage 1, the primary stellar photons 
are followed as they travel from the stellar surface out- 
wards into the surface layers of the flaring disk. In order 
to avoid having to solve a 2-D problem, this is done us- 
ing a "grazing incident angle recipe" (also called "flaring 
angle recipe") which splits the problem into a 1-D verti- 
cal transfer problem at each radius. Photons are inserted 
into the disk at a flaring angle, which is computed self- 
consistently from the disk geometry. As the photons enter 
the surface layers, they get absorbed by the dust and leave 
their energy behind. In the second stage, the re-emission 
and subsequent radiative diffusion of this energy in the in- 
frared regime is modeled by solving a 1-D radiative trans- 
fer problem using the powerful variable Eddington factor 
technique. This technique allows us to solve the transfer 
problem at any optical depth, no matter how high. This 
is a major advantage over the often used (accelerated) 
lambda iteration techniques. It is assumed that the disk 
is passive, i.e. there is no internal production of heat by 
viscosity or any other mechanism. Therefore the above 
procedure suffices to determine the dust temperature at 
every height above the midplane. From this, the dust+gas 
density structure can be determined by integrating the 
equation of vertical hydrostatic equilibrium for a given 
surface density S, under the assumption that the gas and 
the dust temperatures are the same and the gas-to-dust 
ratio is constant (we take it 100). The entire procedure 
is iterated in order to get a self-consistent vertical tem- 
perature and density structure. It typically requires 4 to 6 



iterations for a convergence to within 1%. For more details 
on the model see DZN02. 

In order to include scattering into this model, we re- 
place stage 1 with a Monte Carlo code capable of mod- 
eling isotropic as well as small-angle scattering. In this 
way scattering of the primary stellar photons is included 
in a self-consistent way. The possible scattering of the re- 
cmitted infrared photons can not be treated in this way, 
and will be ignored in the following. In practice, our code 
(and the results presented here) apply to disk atmospheres 
where the grains are relatively small, and their albedo is 
negligible at the wavelengths of the reprocessed radiation 
(i.e., in the infrared). 

The flaring angle a, under which the photons first enter 
the disk's surface layers before they scatter, is given by: 



a = 0.4- 



.R* 
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(1) 



(see e.g. Chiang & Goldreich Il997fk where H S (R) is the 
surface height of the disk, defined as the vertical height Z 
at which the radial optical depth towards the star at stellar 
wavelength is unity. This flaring angle is recomputed after 
each iteration, so that in the end of the iteration procedure 
one obtains a self-consistent solution with the prope r flar- 
ing angle. As was described by Chiang et al. I|200lf) . one 
has to employ a special method to compute the flaring 
angle such that numerical instabilities do not appear. 

2.1. Monte Carlo method 

In the Monte Carlo radiative transfer we let primary stel- 
lar photons impinge on the disk at an angle a with respect 
to the surface. We use a vertical 1-D plane-parallel geom- 
etry, with a vertical spatial coordinate Z ranging from 
Z = (the equator) to Z = Z max . In this approach, the 
incoming photons are assumed to have fiQ = — sin(a) in 
the traditional notation used in 1-D radiative transfer. For 
each frequency bin we use N photon packages. Each pho- 
ton package is assumed to carry an energy (per unit disk 
surface area) of 

In the classical Monte Carlo approach each photon pack- 
age is now followed as it scatters through the medium (ran- 
domly changing it's direction fi at each scattering event), 
until it gets absorbed at some point. Absorption is thus 
treated as a discrete event. By counting the number den- 
sity of absorbed photon packages per grid cell, one can 
determine the energy deposition into each cell. This ap- 
proach, although formally valid, has a serious weakness: 
it produces enormous numerical noise at small optical 
depths. 

A better approach is to treat absorption as a contin- 
uous process, peeling off energy from the photon package 
according to the law: 



dE v 
ds 



= -pK a v E v 



(3) 
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where k" is the absorption opacity, and s is the path 
length along the photon trajectory which obeys the re- 
lation ds = dz/fi. The energy lost by the package is de- 
posited in the cell in which it was lost, while the package 
continues to move on until its energy E has dropped be- 
low some prescribed fraction of the original energy E®, or 
until the package has left the disk. Every time a photon 
package passes through a cell, another portion of energy 
is deposited into that cell, until all photon packages have 
been used up. The total energy (per unit surface area of 
the disk) deposited in the cell constitutes a heating rate q 
per unit volume. 

Scattering is still treated in a discrete way, and the en- 
tire package regularly changes direction due to these dis- 
crete scattering events. At the start of the random walk, 
and after each scattering event, a random number £1 uni- 
formly distributed in the domain (0, 1) is chose n usin g 
the random generator ran2 of Numerical Recipes (^992). 
Using this number, the location of the next scattering 
event is computed according to: 

^ncw — -^current log £1 . (4) 

PK S 

At this location, the Henyey-Greenstein scattering func- 
tion is invoked to compute the change in direction. This 
is first done in the photon's comoving frame. The Henyey- 
Greenstein (1941) function is: 

p(co S e') = l — 2 l ~/ (5) 

2 (1 + g 2 — 2gcos6') A i z 

where g is the forward-peakedness parameter which is a 
property of the dust grains, and 0' is the angle between 
the new photon direction and the old one. Using again the 
random number generator one can determine the 9': 

where £2 is again the random number between and 1. 
The other angle (</>') is simply a random number between 
and 2-7r: </>' = 2tt^. We now employ a rotation matrix 
operation to rotate these angles to the real ones 8, tf>, the 
new /j, being /i = cos 9. 

Using the above procedure, one can compute the heat- 
ing q(z) due to the incidence of primary stellar photons. 
This heating rate then serves as the source term for stage 
2: the radiative transfer of thermal emission from the 
disk. The continuous treatment of absorption in the Monte 
Carlo simulation assures that this function is virtually free 
of numerical noise even for a modest number of photon 
packages (say, N — 100). This method (whi ch is similar 
to the method of Van Zadelhoff et al. [2003) works very 
well at all optical depths. We have tested that the method 
reproduces the results of a Monte Carlo code using dis- 
crete absorption events. 





Monte Carlo code 
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Fig. 1. The fraction of light diffusely reflected off a semi- 
infinite slab as a function of albedo of the slab. Compared 
are the analytic results of Chandrasekhar, the results from 
our Monte Carlo code, and the results from the thermal 
reprocessing computed by our code. 

2.2. Comparing to diffuse reflection theory of 
Chandrasekhar 

The above numerical procedure can be tested against the 
analytic so lutions of diffuse scattering (Chandrasekhar 
Il950 71960^ . Using so called H-functions, the problem of 
photons scattering off a semi-infinite slab with a certain 
albedo u> = K scat / (K a bs + ftscat) can be reduced to a single 
numerical integral. For the case of isotropic scattering the 
ff-function of Chandrasekhar is defined in the following 
way: 

-J- = VT^ + % f -^Hfa'W . (7) 

The values of the 77-function can be quickly found 
by iteratively evaluating Eq.J7|) until concergence is 
reached. Ta bulated values can be found in Chandrasekhar 
l|l950 /1960|L Once this iJ-function is found, one can find 
the reflected fraction of incident flux rj as: 

r ] = l -u J I -J^ H ^ )H(n')dv' , (8) 

where /iq is the incident angle of incoming flux, assuming 
that the flux is a parallel beam. In the case of a disk 
irradiated by a star, this assumption is reasonably well 
justified. The reflection fraction 77 is a number between 
and 1, where rj = 1 means that every photons is eventually 
scattered away without having been absorbed, and r\ = 
means that all photons have been absorbed. 

In Fig.^the reflected fraction of incident flux is shown 
as a function of albedo. Plotted over each other are the 
predictions from Chandrasekhar's theory of H-functions, 
the reflected fraction computed using our Monte-Carlo 
code, and one minus the fraction of the flux that is re- 
processed into the infrared according to our code. All fall 
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reasonably well over each other, so we trust our numerical 
results. 

Fig. ^ shows that for small grazing angle and small 
albedo one has rj ~ oj/2, meaning that the infrared lu- 
minosity of the disk is easily estimated to be Lir ~ 
(1 — oj/2)L*. This result can be understood in terms of 
the Chiang & Goldreich picture of a hot surface layer: 
this layer intercepts all stellar radiation, and redirects half 
of it away from the disk and half towards the disk. This 
redirecting proceeds partly through scattering (u>), partly 
through absorption (1 — uS). This means that a fraction of 
u>/2 of the incident flux is scattered away from the disk 
surface. The other w/2, which is directed towards the disk 
interior, for low albedo will be almost entirely absorbed, 
so that (1 — ui/2) will be reprocessed in the infrared. For 
larger albedo, the downward directed scattered photons 
might be scattered once more, and escape after all. This 
causes the deviation from the i] ~ uj/2 estimate, as seen 
in Fig. [I] 

2.3. Approximate treatment of isotropic scattering 

For isotropic scattering phase function there is an inter- 
esting approximate alternative to th e full M onte Carlo ap- 
proach described above (Calvet et al. ll99li Strittmatter et 
al ll974l) . Using the moment equations of radiative trans- 
fer in the Eddington approximation one can derive an an- 
alytic formula for the mean intensity of the direct + scat- 
tered (but not the re-emitted) photon field J v as a function 
of vertical absorption + scattering optical depth into the 
disk t v . If we define (3 V = -\/3(l — u> v ), and = sin(a), 
where a is the grazing angle of incidence of the imping- 
ing stellar photons, then the mean intensity J v (t u ) can be 
written as: 



J v {t v ) 



Hqu{2 + 3/x )F; 



"47r(l + 2/V3)(l-/M 



47T 



(9) 



4tt(1 - f%$) 

(Calvet et al. ll99lh where F* is the unextincted stellar flux 
at the radius of the disk annulus. Using this expression for 
the mean intensity of scattered + direct stellar radiation 
the heating rate q(z) can be determined: 



47T 



J v {z)pn v dv 



(10) 



This q(z) then serves as the source term to stage 2 of the 
radiative transfer calculation. 

We have compared the models resulting from these 
equations to those resulting from the more precise Monte 
Carlo method described in Sec l2.ll We found that under 
most circumstances, for isotropic scattering phase func- 
tion, the temperature structure derived from the use of 
Eas. (|9ll0|l and that from the Monte Carlo method differ 
by less than 3%. For isotropic scattering, Eq. © therefore 
seems to be a reliable formula. However, it works only 
when the disk has large optical depth at all wavelength 



where the stellar flux is important. Also, it is important 
to note that the equation for the dust temperature pro- 
posed in Calvet et al. (their Eq.14) is not very accurate, 
as it is based on the use of the Rosseland mean opacity. 
As was discussed by DZN02, such an approach may lead 
to a wrong temperature structure. 

3. Single annulus setup 

In order to get an understanding of the effect of scatter- 
ing on the disk structure we start with a single-annulus 
setup. The central star is assumed to have the follow- 
ing parameters: M» = 2M , R« = 2R Q , T.« = 10000K 
(i.e. L» = 36L Q ). These are parameters relevant to Herbig 
Ae stars. For simplicity we assume the spectrum of the 
star to be a blackbody spectrum. The annulus of the flar- 
ing disk we wish to study is located at 1 AU, and has a 
width of 0.01 AU (i.e. it lies between 1.00 and 1.01 AU 
from the central star). The grazing angle (i.e. the angle at 
which the stellar radiation enters the disk's atmosphere) is 
generally computed self-consistently from the disk model, 
but for the simple test case described here we assume it 
to be fixed to a = 0.05. The opacity is take n to b e that 
of a silicate grain of 0.1 fim (Draine & Lee Il984l) . The 
scattering opacity is constructed artificially, in order to 
be able to investigate different values of the albedo. We 
assume that the scattering is zero for A > 2/im, and that 
it is a fixed fraction of the absorption for A < 2/j.m (keep- 
ing the absorption opacity itself unchanged). In this way 
we can freely vary the albedo for the incident stellar ra- 
diation, while the albedo for the reprocessed radiation re- 
mains zero. It should be noted that with this definition 
the total (scattering-l- absorption) opacity for A < 2/im in- 
creases for increasing albedo. The surface density of our 
test setup is defined using the optical depth for the case of 
zero albedo. The annulus has a surface density of S = 0.3, 
and (for w = 0) a visual optical depth in the vertical di- 
rection of Ty = 10. 

For the setup described above we compute a number of 
models varying the albedo and the phase angle. In a first 
set of models we assume that the scattering is isotropic 
and change the albedo in the range 0.5-0.9. In a second set 
of models, we will fix the albedo to be w=0.8 and consider 
forward-peaked scattering by varying the phase angle from 
0.0 to values in the interval 0.5-0.9. In the ISM, typical 
values o f the albedo are of ab out 0.5, with g ^0.5 (Draine 
and Lee ll984t Kim et al. ll994^ . Both parameters, however, 
depend strongly on the grain properties, suc h as si ze and 
chemical composition (e.g. Mishchenko et al. 2000). 

3.1. Varying the albedo 

In Fig.Elthe vertical temperature profile is shown for four 
values of the albedo, from 0.0 to 0.9 and isotropic scat- 
tering. As the albedo increases, an increasing fraction of 
the stellar energy is lost from the heating budget. The ef- 
fect of scattering is to decrease the disk temperature at all 
heights, with the only exception of the very optically thin 
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region on the surface, where the stellar photons reflected 
away from the disk slightly increase the mean intensity of 
the heating radiation. The temperature decrease is larger 
at intermediate optical depth, while the midplanc temper- 
ature drops only slightly, by 6% for w=0.5 and a maximum 
of 23% for w=0.9. It is easy to understand why this effect 
is so small if one loo ks at i t from a two-layer perspective 
(Chiang & Goldreich ll997|) . The disk midplane is heated 
by the emission from the surface layer, part of which is 
thermal dust continuum (1 — oj) and part of it is scat- 
tered light (oj). Only a fraction (oj/2) of the downward 
scattered light has a chance to scatter backwards away 
from the disk, and escape the fate of being absorbed and 
reprocessed. This means that the effect on the midplane 
temperature should be of the order of oj 2 (see W2.2\ . 

Fig. 01 shows the corresponding density profiles. One 
can see the deviation from the gaussian vertical profile 
expected in an isothermal disk. The effect is significant 
at high Z , where the increase of the temperature due to 
stellar irradiation occurs. One can also see that at fixed Z 
the density is lower for higher values of oj. However, the 
effect of scattering on the pressure scale height H p is min- 
imal, since H p depends on the square root of the midplane 
temperature, which changes very little (see Fig- El ■ 

Fig. [3and Fig.OJshow the location of the surface height 
H s , i.e., where the stellar-flux-averaged optical depth (ab- 
sorption + scattering) to the stellar photons entering the 
disk under the angle a is unity. The location of H s , which 
defines the "surface" of the disk, moves upwards as oj in- 
creases, reflecting the increase in total optical depth as 
more scattering opacity is added. Even if for increasing 
albedo the density structure of the disk "shrinks" , this ef- 
fect is counteracted by the increased total opacity, and H s 
in fact moves towards higher Z / R. 

It is interesting to note that intuitively one might ex- 
pect scattering to increase the temperature of the disk 
just below the surface H s , since scattering redirects stel- 
lar radiation deeper into the disk. However, as can be seen 
in Fig. [3 such an increase in temperature does not take 
place. The explanation lies in energy conservation. The 
deep layers of the disk get heated indirectly via emission 
from the surface layer of the disk. Without scattering this 
downward flux is purely thermal emission from the hot 
dust in the surface layer. When scattering is included, the 
total flux that is beamed downward stays the same, but 
now consists partly of thermal emission from the surface 
layer and partly of scattered light at stellar wavelengths. 
The only effect that can be discerned is that, due to the 
different absorption cross sections in these two wavelength 
ranges, the scattered photons are absorbed higher up in 
the disk than the thermal photons. This very subtle effect 
is seen in the oj = 0.9 temperature curve in Fig.|3between 
Z/R = 0.04 and Z/R = 0.08, where the slope changes less 
rapidly than for smaller values of oj. 

Fig. El shows the resulting SEDs for our test annulus 
for the same values of the albedo oj = 0, 0.5, 0.8, 0.9. Only 
the thermal emission from the annulus is shown, not the 
scattered- away radiation which remains at stellar wave- 
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Fig. 2. Vertical temperature profile of the test annulus 
as function of Z/R for different values of the albedo, as 
labelled. The tickmark on each curve indicates the position 
of the surface height H s of the disk, i.e. the height where 
the optical depth (absorption + scattering) with respect 
to stellar photons is unity. 
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Fig. 3. Same as Fig. |2but now for the vertical density 
profile. 



lengths. As one can see, the main effect is to reduce the 
overall infrared flux from the annulus. This is a logical 
consequence of energy conservation: as more radiation is 
scattered away instead of being absorbed-and-re-emitted, 
the total budget of reprocessable radiation is of course 
lower. In this particular example, the fraction of repro- 
cessed (i.e. not scattered) stellar radiation varies from 1 
for oj = to 0.75 for oj = 0.5 and to 0.37 for oj = 0.9. 

At short wavelengths the reduction is stronger than at 
long wavelengths, so that the SED appears redder. At mm 
wavelengths the reduction between oj — and oj = 0.5 is 
about 6% (proportional to the decrease of the midplane 
temperature), while at A = 3/im the reduction is about 
40%. If one quantifies the reddening by the ratio of the 
NIR/IR emitted luminosity, with NIR between 1 and 7 
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Fig. 4. SED of the test annulus for different values of the 
albedo, as labelled. The flux from the star is not included, 
nor is the flux from the radiation that is scattered away 
from the disk. Only the thermal emission is shown. 
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Fig. 6. The effective reduction of the albedo as a result 
of non-isotropic scattering, for different values of uj and 
incident grazing angle a. 
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Fig. 5. The effect of increasing the small- angle scattering 
factor g for constant albedo w on the spectrum of the 
annulus test case. 

/j,m, it reduces from 0.19 for uj = 0, to 0.16 for u> = 0.5 
and 0.09 for uj = 0.9. The reason for this is that the NIR 
emission forms in the upper layers of the disk where scat- 
tering has a stronger effect than near the midplane. At 
very long long wavelengths the disk emission is propor- 
tional to the temperature of the midplane, which, as we 
have discussed, decreases only little. 

3.2. Small angle scattering versus isotropic scattering 

In order to find out which effect a non-isotropic scattering 
phase function has on the structure of the disk and on the 
SED, we take the single-slice setup of the previous section, 
and make runs for different values of g, using the Henyey- 
Greenstein phase function. For a = 0.05 and u> = 0.8 
the differences in the SED are shown in Fig. [3] From this 



figure it can be seen that increasing g has the effect of 
decreasing the influence of scattering on the disk. As the 
g factor increases, the infrared spectrum of the annulus 
also increases, since the scattering becomes less effective 
in scattering radiation away before it can thcrmalizc and 
be reprocessed into the infrared. For increasing g, the SED 
progressively moves towards the SED as computed from 
the case without scattering. 

The reason why this is the case can be understood by 
taking the extreme example of perfectly forwardly peaked 
scattering (g = 1), which is equivalent to no scattering 
at all, since a scattered photon will keep moving in the 
same direction as before the scattering event. For g = 1, 
one can therefore be sure that the results should be iden- 
tical to the case of no scattering, even for high albedo. For 
smaller g scattering becomes important, but the forward 
peakedness of scattering will supress the effect of scatter- 
ing compared to the case of g = 0. This is because the 
forward peakedness of scattering phase function improves 
the chance of still pointing into the disk after the first 
scattering event, but still leaves a chance to get scattered 
away. 

The question now arises: could one simulate the ef- 
fects of small angle scattering by isotropic scattering with 
a reduced "effective" albedo? To answer this we set up 
the following procedure. We start with a calculation of 
the kind described above with a particular u> and g. Then 
we make a new calculation with g = 0, and tune the lu 
until we obtain the same reflection fraction as in the orig- 
inal computation. This 'effective' albedo w e ff can be au- 
tomatically found by placing the model computation as 
a function-call inside a root fin ding r outine (for instance 
zbrent from Numerical Recipes Il992() . We found that for 
any value of the parameters R, a, uj, g we explored, it is 
possible to find a value of uj (called w e ff) that reproduces 
the same SED with g = to within 5%. This effective 
albedo is always uj c g < uj. 
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In Fig. the effective reduction of the albedo 
(i.e. uj e ff/u)) is shown for increasing value of g for different 
values of the incident grazing angle a and original albedo 
ix). One can see that for small u> the reduction is stronger 
than for big uj. This can be explained by the fact that 
high albedo gives a photon more than one chance to scat- 
ter and change its direction. Multiple scattering will erase 
the memory of the original direction where the photon 
came from. One also sees from the figure that a higher 
incident angle a will aggravate the reduction effect. This 
is because for higher a the photon will need a larger scat- 
tering angle in order to be able to escape the disk. In the 
extreme case of a = tt/2 (vertical incident radiation) a 
photon will need to scatter essentially backwards in or- 
der not to get absorbed within the disk. For this case the 
counter effect of g will be the strongest. For the limit- 
ing case of a = 0, small angle scattering will essentially 
change nothing, since even a small scattering angle has 
50% chance of deflecting the photon away from the disk. 

An analytic fitting formula for the effective albedo is: 

uj cS = a,(l - g a ) b (11) 

with 

a = 3.7 uj 2 - 9 - 1.6 a 042 + 2.5 (12) 

b = -0.31w 29 + 1.3a a42 (13) 

This formula has a mean deviation of 10% with the com- 
puted values within the domain of < g < 1, 0.01 < a < 1 
and 0.1 < u < 0.9. 

4. Full disk models 

We turn now to the determination of the effects on the 
SED of a full disk ranging from the inner dust-evaporation 
radius out to a few hundred AU. We have set up a disk 
model for the same Herbig Ae star parameters as in Sec. 01 
mass M* = 2M , effective temperature T* = 10000K and 
luminosity L* = 36L©. The disk surface density is taken 
to obey S(i?) = 400 (R/AU)^ 1 g/cm 2 with an inner radius 
of 1.0 AU and outer radius 300 AU. The mass of this disk 
is A/disk = 0.08A/q. In these models the flaring angle and 
pressure scale height are computed self-consistently, as in 
DZN01. We vary the albedo from uj = to uj = 0.9 and 
assume isotropic scattering. We ignore, for the moment, 
the emission from the inner rim of the disk. 

The resulting SEDs are shown in Fig. As expected 
following the discussion in scattering has the effect of 
reducing the infrared excess over the entire wavelength 
range. The ratio of the reprocessed luminosity to the stel- 
lar one varies from 0.41 for u> = to 0.36 for u> — 0.5 to 
0.19 for uj = 0.9. Roughly speaking, the effect of scatter- 
ing on the infrared excess is to reduce it by a factor of the 
order of (l-o;/2). This is not exact (see Subsection I2.2|) . 
but gives a zero order of magnitude of the effect. 

As in the case of the single annulus setup the reduc- 
tion is stronger at short wavelengths than at long wave- 
lengths. This is more clearly seen in the lower panel of 



Fig. [7| in which the reduction factor is plotted as a func- 
tion of wavelengths. The reduction is practically negligible 
at millimeter wavelengths, while it is as large as 80% (for 
w=0.9) in the near infrared. It is interesting to note that 
the reduction factor is a rather smooth function of wave- 
length, and that it shows only a weak variation around 
the location of the 10 [xm silicate feature. The degree by 
which one sees this effect varies from model to model. It 
is a subtle effect of the precise vertical profile of the tem- 
perature and the density in the disk, and is also linked to 
the slight reduction of the disk's interior temperature for 
increasing albedo. If we enforce the midplane temperature 
of the disk to be unaffected by the scattering (by fixing 
the temperature at z — to the value it had without scat- 
tering), then the 10 /im feature turns up quite strongly in 
the reduction factor. 

Fig.|Hlshows the SEDs of disks irradiated by a T Tauri 
star (TTS). The parameters of the star are taken to be: 
T* = 4000K, R„ = 2R Q , M* = 0.5M Q . The disk has an 
outer radius of 300 AU, an inner radius of i?; n = 0.1 AU, 
a mass of A/disk = 0.008 M©, and a powerlaw index of the 
surface density distribution of —1. The results are qualita- 
tively very similar to the Herbig Ae star case. For uj = 0.5, 
however, an interesting phenomenon occurs. Longwards of 
100 /im the reduction factor for the SED exceeds unity, 
i.e. an inflation instead of reduction. This can be under- 
stood as a result of the increasing height H s of the disk, 
as can be seen in Fig. 1101 The grazing angle a increases, 
enhancing the irradiation. This effect however is rather 
small, and on the whole the infrared excess is reduced by 
the inclusion of scattering. 

In Fig. Elthe same SEDs is shown as in Fig. El but this 
time the scattered light from the disk is included in the 
spectrum (the direct stellar light is still not included). As 
one can see: the flux of the scattered light increases for 
increasing albedo. This is a combined effect of the disk 
becoming more refractory and the increasing height H s . 

In addition to the surface height H s , Fig. ^] also 
shows the pressure scale height at the equatorial plane 
Hp. Contrary to H Sl the pressure scale height varies only 
little as a function of albedo. 

4.1. Inner rim 

For Herbig Ae/Be stars there is evidence that, in comput- 
ing the SED, one should include the effect of the stellar 
radiation impinging frontally on th e disk at the dust sub- 
limation radius (Natta et al. 1200 it Dullemond, Dominik 
& Natta l200ll henceforth DDN01). Since scattering of the 
stellar light in this inner rim may be important, we have 
computed the SED of rim+disk models for various values 
of the albedo. The rim is treated in a way similar to that 
of DDN01. Knowing that the rim is located at the dust 
evaporation radius, we know its temperature (about 1500 
K) and we can compute the pressure scale height of the 
rim. The surface height of the rim (the z above the mid- 
plane at which the rim becomes transparent to starlight) 
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Fig. 7. The effect of increasing albedo on the total SED 
of a passive flaring irradiated disk around a Herbig Ae 
star. The SED is computed at an inclination of i = 45°. 
No inner rim is included. Only the thermal emission from 
the disk is shown. 
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Fig. 8. Same as Fig. for a TTS. 
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Fig. 9. Same as Fig. 00 but now with the scattered light 
from the disk included. The direct light from the star is 
not included. 
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Fig. 10. Run of the pressure scale height H p and of the 
surface scale height H s as function of R for different values 
of the albedo, as labelled. 



is then easily computed. As in DDN01, this height is con- 
siderably higher than the nominal height of the flaring 
disk at that radius, so this rim will in fact cast a shadow 
over the flaring disk. The SED from this disk consists then 
of a component in the near infrared originating from the 
hot inner rim, and the emission in mid- and far-infrared 
originating from the non-shadowed part of the flaring disk 
behind the rim. Scattering is included in the treatment of 
the inner rim using the same 1-D plan-parallel radiative 
transfer techniques used for the disk annuli, but this time 
in the horizontal direction. 

The results are shown in Fig. El The SED is very 
similar to that of models with no rim at long wavelengths 
(in this example approximately for A > 15 (an), where 
the emission comes from regions of the disk which are out 
of the rim shade. The region shorter of about 7 /im is 
deeply affected by the presence of the rim, which repro- 
cesses in this region the intercepted stellar radiation. The 
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Fig. 11. Same of Fig.0 but with the inner rim included. 
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Fig. 12. The effect of varying g on the SED of a passive 
flaring irradiated disk around a Herbig Ae star. In all cases 
it is uj = 0.5. These models include the inner rim. Note 
that the vertical scales are different from previous figures, 
in order to show the effects more clearly. 

effect on the inner rim emission of increasing u> can be seen 
clearly in Fig. 1111 lower panel, and it is qualitatively simi- 
lar to what we have discussed in §3: the emission decreases 
at all wavelengths, with the exception of the very short 
ones, where it is fixed by the requirement that T=1500 K. 
Globally, however, the overall shape of the SED depends 
on u> less than that of models with no rim and, in par- 
ticular, the near-IR excess (A < 7/im) is always roughly 
25% of the total. The total infrared excess of these models 
is always higher than in models with no rim. However, it 
decreases in a similar way with u>, roughly as (1 — lu/2) 
(from 0.61 to 0.50 for u = 0.5 and to 0.32 for u = 0.9). 

We also made a set of rim+disk models for fixed 
albedo (w=0.5), but varying g between isotropic scatter- 
ing and strongly forward-peaked scattering (see Fig. 112(1 . 
Increasing the value of g decreases the effects of scattering 
on the SED, but it does so more effectively for the inner 
rim (near-IR) than for the flaring disk (mid-far-IR). The 
reason is that the inner rim is irradiated under an incident 
angle of i = 90°, while the flaring disk is irradiated under 
a small incident angle. For small incident angles a photon 
is still likely to skim off the surface even for large g, which 
is not the case for large incident angle (see Sec 13. 21 . 

5. Summary and conclusion 

In this paper we investigate the effects of scattering of 
the impinging stellar radiation on the structure and spec- 
tral energy distribution of a passive flaring irradiated cir- 
cumstellar disk. We base our model calculations on the 
detailed vertical structure model of DZN02, in which full- 
fledged radiative transfer is used on a density structure 
in vertical hydrostatic equilibrium. It was shown in that 
paper that the inclusion of full frequency-and-angle de- 
pendent radiative transfer (although in a 1+1-D manner) 



can have profound effects on the vertical structure, most 
notably on the dust temperature at the equatorial plane. 
Here we improve these models by replacing the original 
treatment of irradiation by a Monte Carlo approach that 
can handle small-angle scattering by dust particles 1 

We compared models without scattering to models 
with increasing albedo at the wavelengths of the stellar 
radiation. Note that we keep the absorption cross section 
fixed, so that increasing albedo results in increasing the 
total optical depth (scattering +absorption) for the stel- 
lar radiation. As the albedo increases, the disk gets cooler 
at all vertical heights, with the only exception being the 
very upper layers. This reduction of the temperature takes 
place even at the very optically thick midplane, although 
the effect is smaller than at intermediate optical depth. As 
a result, the disk becomes slightly flatter (i.e. it "shrinks" 
in vertical direction). Still, the surface scale height H s , 
which defines the surface where the stellar light is inter- 
cepted, increases, since the increase of the total opacity 
more than compensates for the physical disk shrinking. 
Disks seen in scattered light will appear thicker than pre- 
dicted by models where scattering is ignored. 

As far as the SED is concerned, increasing the albedo 
has the effect of reducing the thermal emission from the 
disk over the entire infrared wavelength domain. This is 
because scattering reflects part of the stellar light away 
from the disk before it had a change to get absorbed and 
reprocessed into the infrared. To zero order, for u> < 0.5 
the reduction is approximately (1 — w/2), and somewhat 
larger for larger values of the albedo. However, we find 
that the spectrum is more strongly affected at short than 
at long wavelengths, so that the overall SED looks "red- 
der". For a typical case with uj = 0.5 the reduction in the 
near and mid-IR could be as large as 40%, while at mm 
wavelength the reduction is of few percent at most. Models 
computed with different values of the phase function show 
that the effect of scattering decreases strongly when it is 
more forward-peaked. 

In summary, we conclude that the effect of scattering 
of the stellar radiation by grains in the disk has to be con- 
sidered carefully in disk models, unless the albedo is low 
(u> <C 0.5), the scattering is extremely forward-peaked, 
or the required accuracy is not very high. Including scat- 
tering self-consistently in the determination of the disk's 
structure is particularly important when one is interested 
in predicting the disk shape that one sees in scattered 
light, and/or the near and mid- infrared emission. Only at 
very long wavelengths the disk properties (both integrated 
flux and intensity profile) are practically unaffected by the 
fraction of stellar light that is scattered rather than ab- 
sorbed by the disk. 

Acknowledgements. We thank the anonymous referee for many 
useful comments and suggestions which have improved the 
paper considerably. We wish to thank E. Kriigel for kindly 

1 All the modeis presented in this paper can be downloaded 
in numericai form from http://www.mpa-garching.mpg.de/ 
PUBLICATIONS/DATA/radtrans/diskscat/. 



10 



Dullemond & Natta: Scattering in protoplanetary disks 



lending us his 1-D plane-parallel radiative transfer code which Wood, K., Kcnyon, S. J., Whitney, B., & Turnbull, M. 
we used to test our own 1+1D disk structure code. We also 1998, ApJ, 497, 404 

thank G.J. van Zadelhoff for discussions that proved to be use- Wood, K., Wolff, M. J., Bjorkman, J. E., & Whitney, B. 
ful during the development of our Monte Carlo code. CPD 2002 ApJ 564 887 
acknowledges support from the European Commission un- 
der TMR grant ERBFMRX-CT98-0195 ('Accretion onto black 
holes, compact objects and prototars'). 



References 

Calvet, N., Magris, G. C, Patino, A., & D'Alcssio, P. 1992, 
Rcvista Mexicana de Astronomia y Astrofisica, 24, 27+ 

Calvet, N., Patino, A., Magris, G. C., & D'Alessio, P. 1991, 
ApJ, 380, 617 

Chandrasekhar, S. 1950/1960, Radiative Transfer (New 
York: Dover) 

Chiang, E. I. & Goldreich, P. 1997, ApJ, 490, 368+ 
Chiang, E. I., Joung, M. K., Creech-Eakman, M. J., Qi, 

C, Kessler, J. E., Blake, G. A., & van Dishoeck, E. F. 

2001, ApJ, 547, 1077 
Cotcra, A. S., Whitney, B. A., Young, E., Wolff, M. J., 

Wood, K., Povich, M., Schneider, G., Ricke, M., & 

Thompson, R. 2001, ApJ, 556, 958 
D'Alessio, P., Calvet, N., & Hartmann, L. 2001, ApJ, 553, 

321 

D'Alessio, P., Calvet, N., Hartmann, L., Lizano, S., & 

Canto, J. 1999, ApJ, 527, 893 
D'Alessio, P., Canto, J., Calvet, N., & Lizano, S. 1998, 

ApJ, 500, 411+ 
Dominik, C, Dullemond, C. P., Waters, L. B. F. M., & 

Walch, S. 2003, A& A, 398, 607 
Drainc, B. T. & Lee, H. M. 1984, ApJ, 285, 89 
Dullemond, C. P., Dominik, C, & Natta, A. 2001, ApJ, 

560, 957 

Dullemond, C. P., van Zadclhoff, G. J., & Natta, A. 2002, 

A& A, 389, 464 
Kenyon, S. J. & Hartmann, L. 1987, ApJ, 323, 714 
Kim, S., Martin, P. G., & Hendry, P. D. 1994, ApJ, 422, 

164 

Malbet, F. & Bertout, C. 1991, ApJ, 383, 814 
Malbct, F., Lachaume, R., & Monin, J.-L. 2001, A& A, 
379, 515 

Mishchenko, M., Hovernicr, J., & Travis, L., eds. 2000, 
Light scattering by non-spherical particles (Academic 
Press) 

Natta, A., Prusti, T., Neri, R., Wooden, D., & Grinin, 

V. P. 2001, A& A, 371, 186 
Press, W., Teukolski, S., Vetterling, W., & Flannery, B. 

1992, Numerical Recipes in Fortran, Second Edition 

(Cambridge University Press) 
Strittmatter, P. A. 1974, A& A, 32, 7 
Testi, L., Natta, A., Shepherd, D. S., & Wilner, D. J. 2001, 

ApJ, 554, 1087 
— . 2003, submitted 

van Zadclhoff, G.-J. Aikawa, Y., Hogcrhcijde, M. R., & 

van Dishoeck, E. 2003, A& A, 397, 789 
Wolf, S., Padgett, D., & Stapelfeldt, K. 2003, as- 

troph/0301335 



